Realizing Thermalization Via an Exact Simulation of a Finite Quantum Bath 
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We demonstrate the thermalization of a small but arbitrary quantum system by exactly simulating 
the unitary evolution of the system and a much larger but finite system (a "bath") to which it is 
coupled. While the universe, consisting of the system and bath, remains in a pure state, and 
undergoes a reversible evolution, the system equilibrates to the Boltzmann distribution. In this 
general canonical setting, we further show that the underlying mechanism for the relaxation is that 
of "eigenstate-thermalization'. 
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In the last few years considerable progress has been 
made in understanding how thermodynamical behavior 
emerges from quantum mechanics. This fundamental 
question has been approached from a number of different 
directions. Recent results have shown that if a quan- 
tum system is large enough, almost all states will be 
"typical" , and have properties that are close to those of 
micro-canonical averages [TJ [5]. Naturally, such micro- 
canonical averages become the canonical (Boltzmann) 
distribution when the universe is decomposed into a sys- 
tem and bath. Further, it has been shown rigorously, 
and under fairly weak conditions, that large systems will 
induce effectively irreversible behavior in small subsys- 
tems [3-5 , and very recently that system-bath inter- 
actions with random phases will induce thermalization 
in small systems [6]. From another direction, it has 
recently been demonstrated via an explicit simulation 
that a many-body system will thcrmalize to the micro- 
canonical ensemble [7j. (That is, the values of single- 
particle observables are those predicted by the micro- 
canonical ensemble.) This work also revealed that in this 
many-body system, the thermalization was a result of the 
fact that almost every energy eigenstate of the full sys- 
tem gives the micro-canonical values for single-particle 
observables. That is, the thermalization happens at the 
level of the individual eigenstates. This phenomena was 
termed "eigenstate-thermalization" when it was conjec- 
tured by Srednicki in 1994 jS], and also has its origins 
in the work of Shnirelman [S] and Deutsch [10] . Various 
other authors have contributed to this line of work jTTT - 
[14] . and it is further connected to Berry's conjecture [15] . 
In fact, the recent work on typicality by Popescu et al. p] 
and Goldstein et al. [2] also makes eigenstate thermaliza- 
tion plausible — all thing being equal, one would expect 
the eigenstates of a system to be typical states. 

While the numerical results of Rigol et al. [7] demon- 
strated micro-canonical thermalization for a many-body 
system, they did not give canonical thermalization. For 
another many-body example, canonical thermalization 
was indicated to some extent in the simulation by 
Skr0vseth [16], although the results varied significantly 
with the initial state of the system. Our goal here is to 



demonstrate the canonical (Boltzmann) thermalization 
of an arbitrary quantum system by a fixed large quantum 
system (a quantum thermal bath) to which it is coupled. 
Note that in doing so, the temperature must be fixed 
by the properties of the bath alone. We are motivated 
not only by the question of whether this can be achieved 
with feasibly-sized baths, but because the ability to do so 
would open a new avenue for exploring the origin of ther- 
mal behavior in quantum systems. In particular, in our 
work here we are able to show that eigenstate thermaliza- 
tion plays a central role in this very generic setting. We 
also expect that these methods will provide a feasible 
way to explore the dynamics of strongly damped non- 
linear systems, for which approximate master equation 
techniques, the usual workhorse of open-systems theory, 
break down [17] ITS]. 

The Hamiltonian of our universe is given by 

ffunv = tfsys + HgX^ ® V bath + ff bath , (1) 

where H sys is the system Hamiltonian, -ffbath is the bath 
Hamiltonian, X sys is the system coupling operator, y bath 
is the bath coupling operator, and g is a constant setting 
the overall size of the coupling. In what follows we always 
work in the joint energy-eigenbasis of the system and 
bath, so that H sys and -Hbath are diagonal. 

Since the bath must thcrmalize any system, it is the 
properties of the bath, along with those of Y hath , that are 
the key to obtaining thermal behavior. The properties of 
our bath are as follows: 

1) The bath must be chosen to have a density of en- 
ergy eigenstates that increases exponentially with energy. 
This condition is essentially just the usual equilibrium 
thermodynamic assumption — the Boltzmann distribu- 
tion for a small system in contact with a bath results di- 
rectly from the assumptions that 1) the density of states 
of the bath is exponential as a function of energy, 2) that 
the energy of the universe is conserved, and 3) that all 
states of the universe are equally likely. The temperature 
of the bath is given by T = 1/(&b/3), where the energy- 
density of states is D{E) cx exp((3E). By definition, the 
temperature of a thermal bath should not change as en- 
ergy is added (the bath is "big"), which means merely 
that /3 is a constant, independent of E. 
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2) The bath must have a large enough energy range 
that in interacting with the bath, the system can make 
energy-conserving transitions between its highest and 
lowest energy states. For the simulations below, the 
energy range of our system is 3.5/i/i, where /x sets the 
overall energy scale of the simulation. We choose our 
bath to have 5000 states covering the energy interval 
h(j,[3, 20.0465], with /3 = 0.4. 

3) The initial state of the bath must contain a large 
number of energy eigenstates of the universe, which 
translates to a large number of energy eigenstates of the 
bath. We choose the initial state of the bath to be an 
equal superposition, with random phases, of a set of 350 
contiguous bath energy eigenstates spanning the interval 
7yi[12.4, 14.1]. The initial bath state must be chosen so 
that all the states of the system are accessible from any 
initial system-state by exchanging energy with the bath. 
The appropriate window is most easily selected by exam- 
ining the eigenstates of the universe Hamiltonian, which 
we examine further below. 

4) The bath interaction operator, y bath requires some 
complexity — that is, its elements should, at least lo- 
cally, vary in a more-or-less random fashion. For exam- 
ple, interaction operators resulting from two-body inter- 
actions in a many-body system are typically complex in 
this manner. This is simplest to achieve by making the 
off-diagonal elements Gaussian random numbers, with 
unit variance. We set the diagonal elements of y bath to 
zero, so as to minimize their effect on the bath spec- 
trum. We also choose the interaction, along with all 
other contributions to the universe Hamiltonian, to be 
real. This reduces the numerical overhead in diagonaliz- 
ing the Hamiltonian. 

5) The typical elements of the interaction operator, 
gX sys ® y batb mus t be significantly larger than the sep- 
aration between adjacent bath energy levels. This is 
crucial, because to generate energy exchange between 
the system and bath, the interaction must be non- 
perturbative as far as the universe is concerned. Of 
course, in order to preserve the energy levels of the sys- 
tem, the interaction gX sys <g> y bath must be perturbative 
from the point of view of the system. Note that due to 
the exponential spectrum of the bath, the bath states 
at the lower end of the spectrum are necessarily more 
sparse in energy. To minimize the effect of this sparsity, 
we increase the overall size of the elements of the bath 
interaction operator for the low-energy states. While this 
is not strictly necessary, it keeps the interaction non- 
perturbative for the low-energy part of the bath spec- 
trum, and thus improves thermalization for low-energy 
initial states. The precise form of the bath operator we 
use is as follows. If we denote the bath energy eigenval- 
ues as £j (energy increasing with the index i), then the 
matrix elements y^ ath are given by 
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where for i < j, Asi — — £ b , Asj = Ej — Sj-x, 
and Wij is a zero mean Gaussian random number with 
unit variance. Here Asq is the energy gap between the 
lowest two energy states of the bath, chosen to be Ae = 
fiexp(—j3Eo), with E n the bath ground state. For our 
simulations we set g = 5x 10~ 3 /i, and rj = 100 Aeo- 

Now we have determined the structure of the bath, 
it is time to couple it to a small system. To limit the 
numerical overhead, we use a system with just four states. 
Naturally, the bath is required to thermalize any system, 
including any system interaction operator X sys , with the 
only condition that X sys be sufficiently non-commuting 
with H sys . We therefore make an arbitrary choice for the 
system energy levels and interaction operator X sys . The 
energy levels of our the system to be H(i[0.5, 1.5,2.2,4]. 
Denoting the matrix elements of X sys by Xij, we choose 
xi2 = -0.7, Xi3 = 0.3, xu = -0.9, £ 2 3 = —1.2, X 2 4 = 
—0.4 = — X34. The diagonal elements of X sys are set to 
zero, since there is no sense in unnecessarily perturbing 
the system. 

To demonstrate thermalization we must evolve the sys- 
tem for an arbitrarily long time. Obtaining an essentially 
exact evolution for long times can be achieved by per- 
forming a full diagonalization of the Hamiltonian of the 
universe, H unv . Since the Hamiltonian is a (real) 20,000 
dimensional matrix, this diagonalization does require a 
very large RAM. Nevertheless, with currently available 
computing resources, and absolute addressing, this is now 
quite feasible. In fact, we have already diagonalized real 
Hamiltonians that are twice this size, and even larger 
problems are clearly feasible. 

In Fig. [T] we present the results of the simulation, for 
two initial states of the system, being respectively the 
lowest and highest energy levels. Both initial states relax 
as desired to the thermal Boltzmann distribution and re- 
main there, albeit with small fluctuations. Interestingly, 
when the system starts in its ground state, the residual 
fluctuations are larger than when it starts in its highest 
energy state. We will return to this phenomena below, 
which is due to the finite size of the bath. 

We now turn to the question of the eigenstate ther- 
malization hypothesis [?]• If eigenstate thermalization 
occurs, then for each eigenstate of the universe, the re- 
duced state of the system (that is, traced over the bath) 
will be the thermal Boltzmann state. We will refer to the 
population for a system energy eigenstate that results 
from the universe being in a single energy eigenstate, 
as the "eigenstate-value" for that population. These 
"eigenstate-values" are shown in the four plots in Fig. [2] 
In each of the four plots, the horizontal solid line gives the 
Boltzmann population for the respective state. The dark 
(noisy) line gives the eigenstate-values for the population 
as a function of the energy of the eigenstates. While we 
see that the eigenstate-values are in fact not equal to the 
thermal value, since they fluctuated significantly, what is 
striking is that for a broad range of energy in the middle 
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FIG. 1. (Color online) The evolution of the populations of the energy-eigenstates of a nonlinear 4-level system coupled to a 
bath with 5,000 states, for two initial states of the system. The horizontal lines give the populations for the Boltzmann thermal 
distribution at the relevant temperature. In plot (a) the initial state is the one with lowest energy, and in (b) it is the state 
with the highest energy. The insets are expanded versions of the plots for early times, showing the initial relaxation to the 
thermal state. 



of the universe spectrum, the eigenstate values fluctu- 
ate randomly about a mean that is equal to the thermal 
value. This is revealed by taking a moving average of 
the eigenstate-values, and we see that this matches the 
thermal value over a broad range. As the interaction 
strength is increased (or alternatively the density of bath 
states is increased), the fluctuations (deviations) of the 
eigenstate values from the thermal value decrease. This 
can be seen in Fig. [2] from the fact that that fluctua- 
tions of the eigenstate-values decrease as the eigenstate- 
energy, and thus the density-of-states, increases. Thus 
as the size of the bath is increased, the eigenstate-values 
will tend to the thermal values, with the result being 
eigenstate-thermalization. Note that it is quite natural 
that eigenstate-thermalization happens only for a win- 
dow of energy in the center of the spectrum, and not at 
the edges. The reason for this is that if the initial state 
of the bath is on the upper (respectively lower) extreme 
of its energy range, then the system cannot thermalize 
because it cannot transfer energy to (respectively from) 
the bath. The system can only thermalize if the initial 
state of the bath, and therefore the state of the universe, 
is not at the edges of its spectrum. 

Our results also reveal that complete eigenstate ther- 
malization is not required for a bath to succeed in ther- 
malizing an arbitrary quantum system. For our case, the 
fluctuations of the eigenstate-values about the thermal 
values are still large, yet the steady-state fluctuations of 
the system about the thermal values are small. This is 
due to the fact that the initial state of the bath is spread 
over a large number of adjacent energy eigenstates, and 
so the steady-state for the system is given by averag- 
ing over many adjacent universe eigenstates. To achieve 
thermalization, it is thus only necessary that the average 



of a large number of adjacent eigenstate-values correctly 
gives the thermal value, the individual eigenstate-values 
can still have large fluctuations. 

The above analysis reveals why the steady-state fluctu- 
ations of the system are significantly larger when its ini- 
tial state is the ground state (Fig. [1^,) as opposed to the 
highest energy state (Fig. [T]d) . This is because we chose 
the same energy window for the bath in both cases. As a 
result, in the former case, the state of the universe covers 
a lower energy window, a window over which the density- 
of-states is lower, and the fluctuations of the eigenstate- 
values are consequently larger. (The distribution of the 
state of the universe over its energy-eigenstates is shown 
in the dashed box in first plot in Fig. [2] for the two initial 
states of the system). Since steady-state fluctuations of 
the system are due to, and directly proportional to, the 
fluctuations of the eigenstate-values, these will be larger 
when the system starts in the ground state. Further, 
this shows that we reduce these fluctuations if we wish, 
by choosing the initial state of the bath to lie in a higher 
energy window. 

We can also gain insight into why the average of the 
eigenstate-values gives the correct thermal value, inde- 
pendent of the extent to which eigenstate-thermalization 
is realized. The key is that even with the interaction 
turned off (g — 0), the average of the eigenstate-values 
equals the thermal value, precisely because the bath has 
an exponential density of states. All that is required of 
the interaction operator is to preserve this average, and 
a random bath interaction operator achieves this. 

Important open questions to be explored further con- 
cern the class of interaction operators y bath from which 
thermal behavior will emerge, and precisely how it is that 
operators in this class arise from the structure of many- 
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FIG. 2. (Color online) Each of the four plots shows the populations of one of the system energy levels (of which there are four). 
In each plot: The dashed horizontal line is the desired thermal (Boltzmann) value. The dark noisy line is the population given 
by each of the energy eigenstates of the universe, as a function of their energy. The solid light grey line is the moving average 
of the noisy line over a window of 200 adjacent eigenstates. The solid horizontal line is the actual steady-state population when 
the initial state is state 4. The dashed box in the plot for level 1: the distribution of the initial state of the universe, over its 
eigenstates, when the system starts in state 1 (left) and 4 (right). 



body systems. 

Acknowledgments: The authors are indebted to Prof. 
Daniel Steck at the University of Oregon and the Ore- 
gon Center for Optics; this work would not have been 
possible without access to the large memory nodes of his 
parallel cluster, funded by the National Science Founda- 
tion under Project No. PHY-0547926. KJ is supported 
by the National Science Foundation under Project No. 
PHY-0902906. MO is supported by a grant from the 
Office of Naval Research (N00014-09-1-0502). 



[1] S. Popescu, A. J. Short, and A. Winter, Nature Phys. 2, 
754 (2006). 

[2] S. Goldstein, J. L. Lebowitz, R. Tumulka, and N. Zanghi, 

Phys. Rev. Lett. 96, 050403 (2006). 
[3] N. Linden, S. Popescu, A. J. Short, and A. Winter, Phys. 

Rev. E 79, 061103 (2009). 



[4] P. Reimann, Phys. Rev. Lett. 101, 190403 (2008). 

[5] M. Cramer, C. M. Dawson, J. Eisert, and T. J. Osborne, 

Phys. Rev. Lett. 100, 030602 (2008). 
[6] J. Cho and M. S. Kim, Phys. Rev. Lett. 104, 170402 

(2010). 

[7] M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 
(2008). 

[8] M. Srednicki, Phys. Rev. E 50, 888 (1994). 
[9] A. I. Shnirelman, Usp. Mat. Nauk 29, 181 (1974). 
[10] J. M. Deutsch, Phys. Rev. A 43, 2046 (1991). 
[11] M. Feingold and A. Peres, Phys. Rev. A 34, 591 (1986). 
[12] M. Horoi, V. Zelevinsky, and B. A. Brown, Phys. Rev. 

Lett. 74, 5194 (1995). 
[13] V. V. Flambaum and F. M. Izrailev, Phys. Rev. E 55, 
R13 (1997). 

[14] H. Tasaki, Phys. Rev. Lett. 80, 1373 (1998). 
[15] M. V. Berry, J. Phys. A 10, 2083 (1977). 
[16] S. O. Skrovseth, Europhys. Lett. 76, 1179 (2006). 
[17] A. G. Redfield, IBM Journal of Research and Develop- 
ment 1, 19 (1957). 
[18] K. Jacobs, EPL 85, 40002 (2009). 



